function [d] = chebcoeff(t, y, N)
% compute chebyshev coeffs using Axon method.
% Occurs directly from data in x, y
% does n coeffs.
% approximation is of the form f(t) = sum(over i) d(i)*T(i,t)
% note that compared to Press, f(t) = sum(over i) c(i)*T(i,t) -0.5*c(0)
%
d = [];
if(length(t) ~= length(y))
   disp('Chebcoeff - arrays not of same length')
   return;
end

n = length(y); % number of points in the data set to evaluate

for k=0:N-1 % calculate the coefficients
   S=0; R=0; % initialize the sums
   for i=0:N-1 % do the inside sums
      ip = round(i*n/N)+1; % compute the index as sampled
      Tk = chebyshev(k, t(ip), N); % this evaluates Tk(ti)
      R = R + Tk^2;  % sum up the R values
      S = S+Tk*y(ip); % sum up the top sum
   end
   d(k+1) = S/R; % get the nth coefficient
end
return;


function [Tk] = chebyshev(K, t, N1)
% generate recursively the polynomials for chebyshev polynomial order K,
% evaluated at t, given N number of data points
Tk = [];
T = [];
N = N1;
for k = 0:K
   kk = k+1;
   
   switch k   
   case 0
      T(kk) = 1;
   case 1
      T(kk) = 1 - 2*t/(N-1);
   case 2
      T(kk) = 1 - (6*t/(N-2)) + (6*(t^2)/((N-1)*(N-2)));
   otherwise
      T(kk) = (2*k-1)*(N-1-2*t)*T(kk-1) - (k-1)*(N-1+k)*T(kk-2);
      T(kk) = T(kk)/(k*(N-k));
   end
end
Tk = T(K+1);
return;
